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Scattered Data Interpolation 
with Multilevel B-Splines 

Seungyong Lee, George Wolberg. and Sung Yong Shin 

Abstract — This paper describes a fast algorithm for scattered data interpolation and approximation. Multtlevet B-splines are 
introduced to compute a C^-continuous surface through a set of irregularty spaced points. The atgorithm makes use of a coarse-to- 
fine hierarchy of control lattices to generate a sequence of bicubic B-splir>e functions whose sum approaches the desired 
interpolation function. Large F>erformanoe gains are realized by using B-spline refinement to reduce the sum of these furK:tions into 
one equivalent B-spline function. Experimental results demonstiBte that high-fidetity reconstruction is possible from a selected set 
of sparse and irregular samples. 

Index Terms — Scattered data interpolation, nrruttilevei B-splines, data approximation. 



1 iNTRODUCTlbN 

SCATTERED data interpolation refers to the problem of 
fitting a smooth surface through a scattered, or nonuni- 
form» distribution of data samples. This subject is of practi- 
cal importance in many science and engineering fields, 
where data is often measured or generated at sparse and 
irregular positions. The goal of interpolation is to recon- 
struct an underlying function (e.g., surface) that may be 
evaluated at any desired set of positions. This serves to 
smoothly propagate the biformation associated with the 
scattered data onto all positions in the domain. 

There are three principal sources of scattered data: 
measured values of physical quantities, experimental re- 
sults, and computational values [22]. They are fotmd in 
diverse scientific and engineering applications. For exam- 
ple, nonuniform measurements of physical quantities are 
collected in geology, meteorology, oceanography, cartogra- 
phy, and mining; scattered experimental data is produced 
in chemistry, physics, and engineering; and nonunlformly 
spaced computational vedues arise in the output from finite 
element solutions of partial differential equations, and vari- 
ous applications in computer graphics and computer vision. 

These fields require scattered data interpolation to de- 
termine values at arbitrary positions, not Just those at which 
the data is available. This facilitates many useful operations 
for visualizing scattered multidimensional data. For in- 
stance, in medical imaging, scattered data interpolation is 
essential to construct a closed surface from CT or MRI im- 
ages of human organs. In geological applications, the de- 
rived interpolation function facilitates a contour map to be 
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plotted. Computer vision utilizes scattered data interpjola- 
tion to F>erform visual surface reconstruction on sparse 
measurements obtained from feature-based stereo or mo- 
tion. In Image morphing, scattered data interpolation is 
useful for deriving a smooth mapping function from the 
correspondence of feature points between a pair of images. 
That same process is used to achieve Image registration for 
remote sensing applications. 

Despite a flurry of activity in this area, scattered data in- 
terpolation remains a difficult and computationally expen- 
sive problem. The vast literature devoted to ttUs subject 
documents various approaches, many of which suffer from 
limitations in smoothness, time complexity, or allowable 
data distributions [22|. [28]. This paper addresses these 
problems and introduces a very fast algorithm for con- 
structing a -continuous interpolation function from arbi- 
trary scattered data. The algorithm applies an effective B- 
spline approximation technique to a hierarchy of control 
lattices to generate a sequence of functions whose sum 
approaches the desired interpolation function. Large per- 
formance gains are realized by using B-spllne refinement 
to represent the sum of several functions as one B-spline 
function. 

This paper is based on the multilevel B-spline approxi- 
mation technique presented in [33], (34) for image 
morphing. In that work, multilevel B-splines were used to 
propmgate user-specified values at scattered features across 
the image. This paper describes multilevel B-splines in 
terms of scattered data interpolation and significantiy im- 
proves the performance by applying B-spline refinement to 
the control lattice hierarchy. In addition, we consider how 
multilevel B-spline approximation can be used to generate 
an Interpolation function through scattered data points. We 
present a sufficient condition for interpolation and an adap- 
tive representation of the control lattice hierarchy to mini- 
mize memory overhead. 

Although our method applies to multivariate data, we 
limit our presentation to bivariate data for clarity. We as- 
sume that the independent data is 2D and the dependent 
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data is a simple scalar. This is the familiar surface interpola- 
tion problem in which we have data values defined at an 
arbitrary set of positions {(Af^.y^)). To interpolate the data, 
we seek to find a function f[x, ^ which takes on the value 

at (^t.yj, i.e., = r(x^,y^). Together, (x^.y^.zj com- 
prise the scattered bivariate data. The reader should note 
that other names for scattered data are offered in the litera- 
ture, including "random," "nonuniform." "irregular." and 
" arb itrarily located . " 

This paper is organized as follows. Section 2 reviews 
previous work. B-Spline approximation is presented in Sec- 
tion 3. That serves to motivate the discussion of multilevel 
B-Spline approximation, given in Section 4. In Section 5. we 
demonstrate how the approximation algorithm is used to 
p>erform multilevel B-spline interpolation. Section 6 gives 
several examples in various applications. Performance re- 
sults and a comparison to thin plate splines and hierarchical 
B-spline refinement are presented in Section 7. Conclusions 
are given in Section 8. 

2 Previous Work 

There is a vast amount of literature devoted to scattered 
data interpolation. The reader is referred to [22). [28J. [44] . 
[2] for excellent surveys. In this section, we review several 
dominant approaches based on Shepard's method, radial 
basis functions, thin plate splines, and finite element meth- 
ods. We also consider related research in image processing 
and geometric modeling, including multlresolution filtering, 
direct manipulation techniques, and hierarchical B-spline 
refinement. 

One of the earliest algorithms In this field was based on 
inverse distance weighting of data. Developed by meteor- 
ologists and geologists [9]. [10]. it has become known as 
Shepard's method [45]. Shepard defined a C** -continuous 
interpolation function as the weighted average of the data, 
with the weights being inversely proportional to distance. 
This technique suffers from several shortcomings, including 
cusps, comers, and flat spots at the data points, as well as 
undue influence of points which are far away. Furthermore, 
it is a global method reqiiiring all the weights to be recom- 
puted if any data point is added, removed, or modified. 
Franke and Nlelson introduced the modified quadratic 
Shepard's method [21] to address these deficiencies and 
produce <ontlnuous interpolation. 

Another jxjpular approach to scattered data interpola- 
tion is to define the interpolation function as a linear com- 
bination of radially symmetric basis functions, each cen- 
tered at a data point. The unknown coefficients for the basis 
functions are determined by solving a linear system of 
equadons. The coefficient matrix is always fuU. and. for 
large data sets, it may become poorly conditioned and re- 
quire preconditioning [13]. [12]. Popular choices for the 
basis functions include Gaussian, multiquadratlcs [27]. and 
shifted log [19]. [12], [22], [28]. Hardy's multlquadratics are 
among the most successful and applied methods due to its 
simplicity, well-conditioned property, and intuitive results. 
See [19], [20], [41] for numerous tests and comparisons. A 
survey of radial basis functions for image warping can be 
found in [43]. 



Thin plate splines are derived by minimizing the integral 
of the curvatures over the domain among the Interpolation 
functions of the scattered points. They are widely used due 
to their visually pleasing results and stability for large data 
sets. Although they are usually formulated as the solution 
to a variational problem, Duchon has shovm thin plate 
splines to be derived from radial basis functions [11]. The 
numerical solution of thin plate splines can be accelerated 
by the multigrid relaxation technique [41. [5], [47]. An alter- 
native to multigrid relaxation has been presented that uses 
hierarchical basis functions or wavelets to accelerate the 
convergence of an iteradve technique such as conjugate 
gradient descent [46], [25]. Nevertheless, the numerical so- 
lution remains computationally expensive when the inter- 
polation function is computed on a large grid. Recendy. 
thin plate splines have been used to generate smooth warp 
functions for image warping and morphing [35]. [31], [32]. 

Another class of solutions to scattered data Interpolation 
is due to finite element methods. This approach involves 
creating some type of optimal triangulation on the set of 
data points to delimit local neighborhoods over which 
surface patches are defined. These patches are constrained 
to interpolate the original data. There are several criteria 
suggested by Lawson [30] to derive optimal trianguladons 
in which long thin triangles with small angles are 
avoided. Piecewlse linear approximation over the triangu- 
lation is not smooth, achieving only C° -continuity. The 
most common C* method uses the Clough-Tocher trian- 
gular interpolant [7]. [2]. [26]. A related technique was 
proposed in [39]. Trianguladon methods, however, are sen- 
sitive to data distribution, i.e.. long thin triangles cannot 
always be avoided. 

Schumaker [44] proposed a two-stage method that first 
generates a grid of data using any method for scattered 
data interpolation. The second stage applies a standard ten- 
sor product approximation on the grid. The resulting ap- 
proximation may be made to interpolate the original data 
through the use of an iterative process called the delta itera- 
tion, developed by Foley and Nielson [15J. Arge et al. [1] 
proposed <m approximation scheme consisting of three 
stef>s: regularlzation. local approximation, and extrapola- 
tion. They first determine the approximation function by a 
local method at a subset of the ^d points where the data 
density is high. That function is then extrapolated to the 
entire grid by a global method. 

In the image processing community, scattered data in- 
terpolation is necessary to perform reconstruction among 
nonuniform samples. A fine survey of nonuniform recon- 
struction techniques can be found in [23]. A trend in recent 
algorithms has been the use of hierarchical, or multiresolu- 
tion. filtering to extend onto all positions the information 
known only at the sparse and irregular samples. Burt 
proposed hierarchical polynomial fit filtering to yield a 
multiresolution set of low-pass filtered images that can be 
combined to form a smooth surface passing through the 
original data [6]. Mitchell proposed multistage filtering to 
handle highly variable sample density [40]. In that work, 
weighted-average filters are repeatedly applied with ever- 
narrowing low-pass cutoff until the proper bandwidth for 
the display is reached. 
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Recent work in geometric modeling has addressed scat- 
tered data interpoladon. A B-spUne approximation teclinique 
for scattered data was Introduced to direcdy manipulate an 
object modeled by free-form deformation [29]. That tech- 
nique calculates the pseudoinverse of a matrix containing 
B-spline basis function values to minimize the approxima- 
tion error Welch and Witkin propjosed a variational apH 
proach to direcdy manipulate B-spline surfaces with scat- 
tered points or curves (50|. The control lattice is determined 
by minimizing the sum of energy fu actionals to derive a 
smooth interpolating surface. The high computational cost of 
pseudoinverse and energy minimization solutions render 
these methods prohibitively expensive for large data sets. 

Forsey and Bartels proposed hierarchical B-spline re- 
finement for object modeling [16J. They augmented B-spline 
approximation with that technique to interpolate a grid of 
data using a control lattice hierarchy [17], (18). Interpola- 
tion is achieved by successively improving the approxima- 
tion at a coarse level with a correction term from the next 
finer level. Although their method is similar in spirit to the 
multilevel B-spline approximation presented in this paper, 
their method cannot handle scattered data. 

3 B-Spune Approximation 

Recentiy. a B-spline approximation technique has been 
proposed for image morphing [33], [34]. In this section, we 
elaborate on that technique in terms of scattered data inter- 
polation and present the details of the algorithm. 

3.1 Basic Idea - _ 

Let n = ((x.y}| 0 < X < m, 0 < y < n} be a rectangular do- 
main in the xy-plane. Consider a set of scattered points 
p= [{x^,y^,z^)} in 3D space, where (x^.yj ^ ^ P^^^ ^• 
To approximate scattered data P. we formulate approxima- 
tion function fas a uniform bicubic B-spline function, which 
is defined by a control lattice O overlaid on domain CI. 
Without loss of generality, we assume that <I> is an (m + 3) x 
(n + 3) lattice which spans the Integer grid in ft (Fig. 1). 
Later, we shall consider the effect of different lattice sizes 
on the approximation function. 

Let ipy be the value of the ij-th control point on lattice 

located at (i.j) for i = -l. 0 m + 1 and J = -l. 0 /i+ 1. 

The approximation function f is defined in terms of these 
control points by 

/(*.j')=ii:^m(o«(,.*H/.o- 

where / = W - hj = I yj - 1. 5 = X - [xj. and t = y-LXJ* ^ 
and Bi are uniform cubic B-spline basis functions defined as 

^(0 = (l-t)V6. 
Bi(t) = (3f^-6f^ + 4)/6. 

Bj(t) = (-3f^ + 3f^ 3f + l)/6 . 

B,(f) = f'/6. 

where 0 < t< I. They serve to weigh the contribution of 
each control point to f{x, y) based on its distance to (x. y). 
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Fig. 1 . The configuration of control lattice 

With this formulation, the problem of deriving function f is 
reduced to solving for the control points in <I> that best ap- 
proximate the scattered data in P. 

To determine the unknown control lattice <D. we first 
consider one data point (x^.y^. ^c) ^ ^ From (1), we know 
that function value f{x^,yc) relates to the sixteen control 
points in the neighborhood of (x^.y^.) • Without loss of gen- 
erality, we may assume that 1 < x^.y^ < 2. Then, control 
points for / = 0. 1. 2, 3. determine the value of fat 
{x^.yc) . For function fto take on the value z/at (x^.yj , the 
control points must satisfy 

3 3 

where = Bj^(s)fi^(0 and s = x^ - 1. r = y^ - I . 

There are many values for the 0jy*s that satisfy (2). We 
choose one in the least-squared sense that minimizes 
5^A=o ^u- ™s tends to minimize the deviation of /"from 
zero over the domain ft. This property will prove useful 
when we apply the B-spllne approximation to a hierar- 
chy of control lattices in Section 4. Simple linear algebra 
using pseudoinverse is used to derive the solution [29]: 

9u ~ 



(3) 



In this solution, the control points near (x^.y^) get larger 
values because they are associated with larger w)y. The re- 
sulting fuiKTtion /has the value at (x^.y^) and tapers off 
smoothly. 

Now, we consider all the data points in P. For each data 
point, (3) can be used to determine the set of 4 x 4 control 
points in its neighborhood. These neighborhoods may 
overlap for sufficlentiy close data points. Fig. 2a depicts two 
such data points, P| and . They may assign different val- 
ues to several shared control points. In general, we resolve 
multiple assignments to a control point ^ by considering the 
data points in its 4 x 4 neighborhood (Fig. 2b). Only these 
points may influence the value of ^ by (3). We call this set 
of points the proximity data set of 0, 

Let Py be the proximity data set of control point such 

that 

Py={{^c^yc^^c)^P\ i-2^x,<i+2.J-2<y,<7+2}. 
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(a) overiapping control points (b) proximity data set 

Fig. 2. Positional relationship between data and control points. 

For each point (x^.y^, z^) in/^^, (3) gives a different value 



0c = 



(4) 



where = w„ = B,{s)Bj(t), A = (i + 1) - [xj. 7 = (/ + 1) - 

[y^ Jr s = - J, f = - [yd ' To compromise among the 
vcilues. <f>y is chosen to minimize error e{(py) = 

2^c(^c0//~ ^c0c)^* The term (w^0jj - ^c^c) ^ dilTerence 
between the real and expected contributions of to furK- 
tlon fat {x^rYc) • other words, it is the approximation er- 
ror due to <l>ij, assuming that the other control points sur- 
rounding (x^^Yc^z^) have their values determined by (3) 
using that data point. By differentiating the error ei^y) with 
respect to 0^. we get 

The proximity data set Py is the set of points in P at 
which control point 0^ has an influence on function / 
When Py contains several points. (5) provides a least- 
squared solution to <l»y which minimizes a local approxima- 
tion error. When Py consists of only one point. (5) reduces 
to (3), leaving no approximation error. When Py is empty, 
however, (py has no influence on /(x^.y^) for any data point 
(x^.y^'-^c) *n ^- This implies that ijty may be assigned any 
arbitrary value, such as 2:ero or the average of the z^*s, with- 
out affecting the approximation error. In that case, we as- 
sign zero to <fty to make function / tend to zero in its neigh- 
borhood. This will prove useful for multilevel B-spline ap- 
proximation in Section 4. 

3.2 Algorithm 

To determine the control lattice <I> from the data points in P 
using (5), It is not necessary to expllcidy identify the proxim- 
ity data set for each control point Since each data point in P 
iiifluences a set of 4 x 4 neighboring control points in <P, it 
belongs only to the proximity data sets of those control 
points. Hence, we can efficiendy accumulate the numerator 
and denominator of (5) for each control point by considering 
each data point in turn. The value of a control point is then 
obtained by division if the denominator is not zero. A null 
denominator occurs only when a control point has an empty 
proximity data set. In that case, we assign zero to the control 



point The following pseudocode outlines this B-spline ap- 
proximation method, which we denote as the BA algorithm. 

BA Algorithm 

Input: scattered data P = {(x^.,y^, z^)} 
Output: control lattice <I> = {<l>y] 
for all ij do let Sy = 0 and cOy = 0 
for each point (x^.y^. 2^-) in Pdo 

let/=[xJ-land7 = LyJ-l 

let s= x,-|_xjand t = y^-ly^j 

compute Wjy's and I^q ZJ^q 

for*. J = 0, 1.2. 3 do 

compute with (3) 

add wi/t>,f to 

add to <o^,,,y^,f^ 

end 

end 

for all ijdo 

if Wy * 0 then 

compute i^y - Sy/COy 

else let <f»y = 0 

end 

The time and space complexity of the BA algorithm is 
0{p + mn). where p is the number of data points, and {m + 3) 
X (n + 3) is the size of the control latdce. Although the con- 
trol p>olnt values are determined locally, we minimize the 
approximation error so that the resulting function properly 
reflects the scattered data. Fig. 3 shows an example. Figs. 3a 
and 3b, respectively, show data P and the B-spllne function 
/derived by the BA algorithm for m = n = 16. Notice that / 
nicely approximates P at densely distributed data points 
and interpolates Pat isolated points. 

The density of control lattice O overlaid on domain U di- 
rectly affects the shape of approximation function / If we 
select lattice O of size (n/ + 3) x (rf + 3) . instead of {/n ■»- 3) x 
(n + 3), then the ijth control point in the new latdce is lo- 
cated at (i-^ . in n. As O becomes coarser, the proxim- 
ity data set of each control p»oint covers a larger number of 
points in P. This causes many data points to be blended 
together to yield a smoother shape for /at the expense of 
approximation accuracy. However, as <P becomes finer, the 
influence of a data point Is limited to smaller neighbor- 
hoods. This enables P to be more closely approximated, 
although / will tend to contain local peaks near the data 
points. These characteristics are evident in Figs. 3c and 3d. 
where m = n = 8 and in = /i = 32, respectively. 

The BA algorithm runs very fast even when the number 
of data points is large. Given 5,000 data pmints. for instance, 
the control lattice for the approximation function can be 
derived in 0.15 second on a Sun SPARC 10. Furthermore, 
since B-splines have local support, only a small neighbor- 
hood in the control lattice needs to be updated when a data 
point is added or deleted. This modification is easily com- 
puted as long as we keep the variables Sy and aiy. 

The approximation function / from the BA algorithm is 
-continuous because it is a bicubic B-spline surface 
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(a) given data (b) m = n = 16 



(c) m = n = 8 

Fig. 3. Examples of B-spline approximation under difFerent control lattice 

generated by control lattice <l>. Since * can be derived very 
quickly from the data points, most of the time for comput- 
ing function /is taken by the evaluation of /from <t> on do- 
main To accelerate this evaluation, we use the forward 
difference method [14). [51] or table lookup for basis func- 
tions depending on the control point spacing. It takes 0.77 
second for a Sun SPARClO to obtain function fon a K024 x 
1,024 grid from die data points in Fig. 3. 

4 Multilevel B-Spune Approximation 

A tradeoff exists between the shape smoothness and accu- 
racy of the approximation function generated by the BA 
algorithm. In this section, we present a multilevel B-spline 
approximation algorithm to circumvent this tradeoff. The 
resulting function simultaneously achieves a smooth shape 
while closely approximating the given data P. The algo- 
rithm makes use of a hierarchy of control latdces to gener- 
ate a sequence of functions ^ whose sum approaches the 
desired approximadon function. In the sequence, a function 
from a coarse lattice provides a rough approximation, 
which is further refined in accuracy by functions derived 
from finer lattices. We further optimize this process by us- 
ing B-spline refinement to reduce the sum of these func- 
tions into one equivalent B-spline function. 

4.1 Basic Algorithm 

Consider a hierarchy of control lattices. <I>0'*^i 

overlaid on domain £2. We assume that the spadng between 
control points for is given and that the sj>acing is halved 
from one latdce to the next. Therefore, if <I>i is an (/n + 3) x 
(n + 3) lattice, the next finer lattice <1>a+i will have {2m + 3) x 
(2n + 3) control points. The position of the ij-th control point 
in Ojt coincides with that of the (2i, 2J)-th control p»oint in <t>jt+i. 

The multilevel B-spline approximation begins by apply- 
ing the BA algoritlun to P with the coarsest control lattice 
4>o. The resulting function ^ serves as a smooth Initial ap- 
proximation that possibly leaves large discrepancies at the 



(d) m = n 32 

resolutions. 

data p>oints in P. In particular. ^ leaves a deviation 

A*Z^ = - t (^c'/c) 1^1^^ Uc'/c* ^c) ^ ^ 

finer control lattice 0| is then used to obtain function f, 
that approximates the difference fj = {(x^.y^, A*z^)}. Then, 
the sum ^ + ^ yields a smaller deviation A^z^ = 
Zc- tUr.yJ- ^Uc.yc) for each point U,.y,,z,) in P. 

In general, for a level Jtln the hierarchy, we derive func- 
tion by using control lattice <I)jk to approximate data = 
{(x,.y,,A*z,)}. where A*z, = z, - Sj:o' ^{^c. ^c) = A*'*^c " 
fj_,(x^.y^). and iSz^ = z^. This process starts from the 
coarsest lattice 4>o continues Increment^ly to the finest 
lattice The final approximation function f Is defined as 
the sum of functions f^, i.e.. f = 4- Note that only the 

coarsest lattice is applied to the original data P to derive 
the global shap>e of function f. All successively finer lattices 
serve to approximate and remove the residual error. In this 
manner, we have an incremental solution for function f 
that yields a smooth and close approximation to P. The 
following pseudocode outlines the basic algorithm for 
multilevel B-spline approximation, which we denote as 
the basic MBA algorithm. Note that a hierarchy of control 
lattices is sufficient to represent function f because each 

can be represented by O*. and f is the sum of the ^'s. 

Basic MBA Algorithm 

Input: scattered data P = {(x^.y^, z^)} 

Output: a control lattice hierarchy <&0' 

let A=0 
while A < /i do 

let P^ ={{x„y,.A*z,)} 

compute <l>^ from by the B A algorithm 

compute A**'zc = A*z^ - f^{x^,y^ for each data point 

letilr=ir+ 1 

end 
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(c) fT^ = ^0 = 16,/n, = ftft - 64 (d) f7^ = fT^ = = = 8 

Fig. 4. Examples of multilevel B-spline approximation under different resolutions for the coarsest and finest control lattices. 



Let p be the number of data points in Pand let (ni + 3) x 
(n + 3) be the size of the finest control lattice <t>fr The num- 
ber of control points in lattice ^i^is a quarter of that in the 
next finer lattice <l>jt+j- Hence, the time complexity of the 
basic MBA algorithm is 0{p+ mn) 0(p+ mn) -^...-i- 
0(p + mn) = 0(hp + 4 mn) . The _ space complexity is 
0{p + -J mn) because we have to store all the control lattices 
in the hierarchy. A function from the basic MBA algorithm 
is -continuous because it is the sum of -continuous 
B-spllne functions. Fig. 4b shows an example. The algo- 
rithm is applied to the data shown in Fig. 4a. which is the 
same as that in Fig. 3a. Notice that multilevel B-spline ap- 
proximation generates a much smoother and more accurate 
function than B-spline approximation. 

In the multilevel B-spline approximation, the density of 
the coarsest lattice Oq determines the area of influence of a 
data point on function f. Larger spacing between control 
points serves to merge the Influences of data points over 
large areas to produce smoother function shapes. On the 
other hand, the density of the finest lattice <I>a controls the 
precision with which f approximates the data points. When 
Oft is sufficiently fine compared to the data distribution, f 
can interpolate the data without an approximation error- 
Fig. 4 shows several examples. In Fig. 4a. the given data is 
the same as that in Fig. 3a. Let (m^ + 3) x (n^ +3) be the size 
of control lattice O^. In Fig. 4b, = r^j = 1 and /lift = = 64 . 
In the figure, the approximation function has a smooth shape 
and interpolates the data points. Fig. 4c depicts the function 
when gets finer = = 16 . In Fig. 4d, a coarser O,, with 
mf, = = S results in an approximation error. 

In practice, we use the same horizontal and vertical 
spacing of control points in a control lattice. The sizes of 
control lattices in the hierarchy are related to the aspect 
ratio of domain O. For the coarsest control lattice Oq. we 
generally choose = 1 or = 1 , depending on the aspect 
ratio, to make Oq as coarse as possible. Then, the basic MBA 



algorithm is applied through the hierarchy until we reach a 
sufficientiy fine control lattice 0ft for which the maximum 
error falls below a specified threshold. This results in an 
approximation function that satisfies our goal for smooth 
function shap>e and accuracy. 

4.2 Optimization with B-Spiine Refinement 

The basic MBA algorithm generates a control lattice hier- 
archy that represents the approximation function f. To 
evaluate f, we must determine function ^ froni control 
lattice Ojt for each level A, and add them over domain Q 
(Fig. 5a). This introduces a significant overhead in compu- 
tation time if / has to be evaluated at a large number of 
points in H, We propose to address this problem by pro- 
gressively applying B-spline refinement to the control lat- 
tice hierarchy- This allows / to be represented by one B- 
spline function rather than the sum of several B-spline func- 
tions. Consequentiy. the computation of is limited to the 
small numba* of control points in rather than all the 
points in Q. 

Let i^O) be the B-spline function generated by control 
lattice 0 and let denote the size of O. With B-spline re- 
finement, we can derive control lattice from the coarsest ^ 
lattice <l>o such that J*\OiJ = and |<I>o| = |<I>,|- Then, the 
sum of functions ^ and ^ can be represented by control 
lattice *f , which results from the addition of each corre- 
sponding f>air of control points in <t>J and <Di. That is, 
n^i) = ^i = 4 + /i.wherey, =0J+<t>,. 

In general, let = Zf^o ^ be the partial sum of functions 
fj up to level ir in the hierarchy. Suppose that function g^,, 
is represented by a control lattice ^P^., such that 
l^'j^ j| = |ct>j^ i| . In the same manner as we computed 
above, we can refine ^'^^i to obtain ^'^^j. and add 4';_, to 

to derive H'* such diat F{H'^) = and = l*^! . That is, 
4*^^ = *f;_i + <I>4. Therefore, from go = ^ and = Oq- we 
can compute a sequence of control lattices to progres- 
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Fig. 5. Approximation function evaluation in the MBA algorithm. 
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sively derive the control lattice H'a for the final approxima- 
tion function / = g^. This process Is depicted in Fig. 5b. 

There are many methods for refining a control lattice 
into another so that they both generate the same B-spllne 
functions [8], [37J. [31, [38], [36]. In this paper, an (m + 3) x 
(n + 3) control lattice 0 is always refined to a (Ztji + 3) x 
(2n + 3) lattice <P' whose control point spacing Is half as 
large as that of <P. With this restriction, we can simplify the 
refinement algorithm for B-spline curves in [38] to derive 
from <t>. Let <j>y and be the ij-th control points In 0 and 
<I>'. respectively. Then, the position of control point 0'2/.2/ 
O' coincides with that of control point in 4>. The vjuues 
for the control points in ^ are obtained finm those in O by 



P21.2J 



+ 6{0,_, J + 0, j_i + + 360;jl, 

= [^i-ij + ^1-1.7+1 + ^if 1.7 ^/+i.7+i 



+ 6(0/.7 + */.7+|)l» 
1 , 

^2i+i.27 = 16 l^i.7-1 + ^/.7+l */+i.7-i ^i+i,7+i 
+ 6(0,7 + 0,^,,^)]. 

^2/+i.27+i = 4" [*/7 ^/J+i ^Mij + *i+L7+i'* 



To generate control lattice ^ji, from data P, we do not 
have to keep the elements of <D4, H'J . and for all k. If 
we apply the B-spllne approximation and refinement to- 
gether at each level, H'/, can be derived by traversing the 
control lattice hierarchy from the coarsest to the finest lev- 
els. In the traversal, only one variable for each of the data 
and control lattice sequences is sufficient to manage the 
computation and intermediate result. This technique, which 
we call the MBA algorithm, is outlined in the following 
pseudocode. In the algorithm, P- f(0) denotes the up- 
dated data {(x,.y^.A**'zJ}. where P = {(x^.y^, A*z^)} and 

MBA Algorithm 

Input: scattered data P = {U^. y^. zj} 
Output control lattice H* 
let 0 be the coarsest control lattice 
let ^' = 0 

while 0 does not exceed the finest control lattice do 
compute 4* from P by the B A algorithm 
compute P = P - fI<D) 
compute *f = 4*' + <& 
let 0 be the next finer control lattice 
refine V into whereby F^i") - fl[^) and = |<ti| 

end 

Let p be the number of data points in P and let (m + 3) x 
(n + 3) be the size of the finest control lattice <Dj^ The 
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number of operations required for B-spline refinement is 
linear in the number of control points. Hence, the time com- 
plexity of the MBA algorithm is 0{hp + mn) . Since we use 
only one variable for each <I>4 , 4^^ . 4*^ . and , the^ space 
complexity of the algorithm is 0(p + mn). If we note that 
the depth of the hierarchy is much less thcin the numlier of 
control points, then the MBA algorithm runs with virtually 
the same time and space complexity as the BA algorithm 
for the finest lattice. However, the MBA algorithm gener- 
ates a much smoother function shape with the same ap- 
proximation error. 

Suppose that the approximation function f must be known 
at M X N grid points on domain 1 i. The processes for evalu- 
ating /in the basic and the optimal MBA algorithms are illus- 
trated in Fig. 5. The basic MBA algorithm requires 0{hMN) 
time to evaluate fon Q since it takes 0{MN} time to evaluate 
each function ^ across all h levels of the hierarchy. The opti- 
mal MBA algorithm, however, requires mn) time to ap- 
ply the B-spline refinement to the control lattice hierarchy 
and 0(MN) time to evaluate the final B-spline function on fi. 
Although both algorithms generate the same function f, the 
optimal MBA algorithm with B-spline refinement realizes 
large computational savings, especially when MN>Tnn. 
When f is evaluated on a 1.024 x 1,024 grid, it takes 0.92 sec- 
ond on a Sun SPARC 10 for the optimal MBA algorithm to 
generate the function shown in Fig. 4b. 

5 Multilevel B-Spune Interpolation 

The MBA algoritlim generates an approximation function f 
that passes near the data points P. but not necessarily 
through them. We now consider how the MBA algorithm 
can be used to interpolate the data. Recall that function 4. 
for level /r > 0 in the hierarchy, is derived to approximate 
and remove the residual error A*z^ at each data point. The 
fined function f is made to interpolate data P once this error 
goes to zero at some level k. In this section, we present a 
sufficient condition for control lattice <1>4 to generate a 
function that removes any residual error. We also 
present an adaptive representation for the control lat- 
tice hierarchy to minimize memory overhead when the 
finest control lattice is required to be very dense. 

5.1 Sufficient Condition for interpolation 

Let Pi =Ui.yi.Zi) and ft ^i^^Yz^h) ^0 points in 
data . Without loss of generality, we assume that has 
the same horizontal and vertical spacing between control 
points. We define the distance between p, and ft as 

njaA^||^i^J-[^^J,|^J-|^^Jj, where s is the control point 

spacing. This distance represents the maximum number of 
horizontal or vertical lattice lines in 4)* that lie between Pj 
and ft after they have been projected onto domain £1. We 
use it to define the interpolation property in terms of the 
control lattice density and the data distribution. 

Let d be the minimum distance among all pairs of data 
points in /J . If d > 4 , then no two data points share a con- 
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Fig. 6. Examples of the interpolation and approximation cases. 

trol point in their 4x4 neighborhoods. In that case, each 
control point in contains at most one data point in its 
proximity data set. When the BA algorithm is applied to 
with Ojt. (5) reduces to (3) for all control points. This results 
in an interpolation function ^ that satisfies (2) at each data 
point in /J . If d < 3, however, there is a control FK)int in 
that has at least two data points in its proximity data set. In 
that case, the efiiects of those points are blended together to 
determine the value of ^ so that function only approxi- 
mates them. 

Figs. 6a and 6b show examples of the interpolation and 
approximation cases, respectively. Note that the interpola- 
tion property depends on the number of control lattice lines 
between data p>oints, not the actual distance between them. 
This is depicted in Fig. 6 where points and P2 have the 
same distance but different interpolation properties. 

5.2 Adaptive Control Lattice Hierarchy 

The interpolation condition given in Section 5.1 demon- 
strates that the MBA algorithm can interpolate data when 
the control point spacing in the finest lattice becomes 
sufficientiy small to satisfy d > 4 . It is possible, however, 
that a single pair of close data points may require O/, to be- 
come very dense even though all other data points are 
sparsely distributed. This can impose large memory over- 
head in computing the interpolation function. To avoid this 
drawback, we prop»ose an adaptive representation for the 
control lattice hierarchy that only stores those control points 
that lie in a 4 X 4 neighborhood about each data point 

When the B-spline approximation is applied in the MBA 
algorithm, control jioints with empty proximity data sets are 
assigned the value zero and do not contribute to the final 
function f. Therefore, to derive f, it is sufficient to maintain 
only those control points that have data points in their 
proximity data sets. We can then represent a control lattice as 
a set of necessary control points rather than a lattice of all 
control points. As successive B-spline approximations pro- 
ceed across a control lattice hierarchy, marry of the control 
points In the finer lattices have empty proximity data sets. 
Hence, the set representation of a control lattice makes it pos- 
sible to save a lot of storage, especially for finer lattices. 

At a level k in the hierarchy, the values of the necessary 
control points in lattice can efficiently be computed by 
modifying the BA algorithm. The variables Sy and eOy are 
allocated in linear arrays instead of two dimensional arrays. 
When data points are considered in sequence, their influ- 
ences on the neighboring control points are stored in the 
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linear arrays, together with the control point indices. The 
arrays are then sorted by index and merged together to 
compute (5) for control points with nonempty proximity 
data sets. The control points not included In the arrays are 
assumed to be zero because they have empty proximity 
data sets. The following pseudocode reflects this modifi- 
cation, which we refer to as the adaptive BA algorithm. 

Adaptive BA Algorithm 

Input: scattered data P = {{x^,yc* z^)) 
Output: a set of control points {0^} 
letr=0 

for each point (x^^Yc, ^c) In Pdo 

set ij, 5. r. Wjy , L L as in the BA algorithm 
for 7=0. 1,2. 3 do 

compute with (3) 

store wf^jy and index (/ + A, J + /) to 5^ 
store vv^ and index (i + Ar. J + /) to 
ler r = r + 1 

end 

end 

sort 5, and by index 
for each index (i, j) in do 

set fj. ^ so that index{o}^) = (ij) for r, < r < 

compute (t»ij = ZH,. ^r/^%r, 

end 

Let p be the number of data points in P. The time com- 
plexity of the modified BA algorithm is 0(p log p) because 
arrays and cu^ : quire sorting. The space complexity is 
0(p). The time and space complexities of the algorithm de- 
pend on the number of data points, not the size of the con- 
trol lattice. Hence, the algorithm will be useful at a finer 
level of the hierarchy where the size of the control lattice is 
much larger than the number of data points. Note that the 
set of stored control points need not form a sublattice be- 
cause a lattice structure is not required for the algorithm to 
handle the control points. 

With the set representation of a control lattice, the final 
function f of the MBA algorithm must be represented by 
the sum of the B-spline functions derived at each level of 
the hierarchy. Function f cannot be converted into one 
B-spline function on the finest control lattice because we 
do not maintain aU of the control points In that lattice. In 
that case, we cannot benefit from the optimization with 
B-spline refinement in evaluating f although we can 
save a lot of storage in deriving f, A reasonable solution 
lies in maintaining the full control lattice at the coarser 
levels so that B-spline refinement can be applied, and then 
switching to the adaptive approach for finer control lat- 
tices. This realizes performance gains with minimal mem- 
ory overhead. 

6 Applications 

Multilevel B-spline interpolation is well suited to any ap- 
plication that can be cast as a surface fitting problem. In this 
section, we demonstrate its use in three diverse applica- 



tions: image reconstruction, image warping, and object 
reconstruction. The scattered data values In these cases con- 
sist of image Intensities (ID), positional constraints (2D), 
and surface points (3D), respectively. These examples dem- 
onstrate the use of the MBA algorithm on multidimensional 
data for various applications. 

6.1 Image Reconstruction from Nonuniform Samples 

Consider a set of nonuniform samples taken from an imaged 
Reconstruction refers to the interpolation necessary to re- 
cover the image from its samples. A robust solution to this 
problem is of great practical importance in image compres- 
sion. Significant data reduction is possible by representing a 
full image with only a few scattered samples. We now ad- 
dress this problem by interpreting a grayscale image as a 
surface, where the value of each pixel represents its height. 
Image reconstruction from nonuniform samples is then cast 
into a surface fitting framework that can be solved with the 
MBA algorithm. 

Consider the image in Fig. 7a. To derive a set of nonuni- 
form samples, we first apply the Sobel operator [24] to the 
image and threshold the result to identify edges. These 
edge contours consist of an important set of. pixels which 
capture the visual details of the image. Due to their rather 
arbitrary definition and distribution, edges alone are not 
sufficient to recover the image. Therefore, we also sample 
the image on a coarse regular grid to prop>erly reconstruct 
the parts of the image where there are no nearby edge 
points. For the purpose of this example, we used a uniform 
sampling rate of one sample per six pixels along each direc- 
tion. The positions of the sampled pixels are shown in Fig. 
7b. The MBA algorithm is applied to the samples to recon- 
struct the surface in Fig. 7c, The reconstructed image that 
correspKjnds to this surface is shown In Fig. 7d. 

Large data reduction is achieved here. The dimensions of 
the original image in Fig. 7a is 360 x 243. The number of 
sample points in Fig. 7b is 7.301. of which 4.984 lie upon the 
edge contours. Therefore, sample points constitute only 8.3 
percent of the total number of pixels in the original image. 
In the MBA algorithm, we used n\, = 2 . = I and 
rUfj - 360 . /ift = 243 fcMT the size of the coarsest and finest 
control lattices, respectively. Fig. 7d demonstrates that the 
MBA algorithm can adequately reconstruct an image from 
a proper set of samples. The original and reconstructed im- 
ages are visually indistinguishable. 

6.2 Image Warping 

Image warping deals with the geometric transformation of 
digital images. It requires a warp function to establish a 
spatial transformation between the input and output im- 
ages. Depending on the application, warp functions may 
take on many different forms. Simple transformations may 
be specified by analytic expressions including afflne. projec- 
tive, bilinear, and polynomial transformations. More so- 
phisticated warp functions that are not conveniendy ex- 
pressed in analytic terms can be determined from a sparse 
collection of feature points which are displaced to define a 
transformation. The displacements for the remaining points 
are evaluated through interpolation. 
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(a) original image (b) sample points 



(c) reconstructed surface (d) reconstructed image . 

Fig. 7. The MBA algorithm is applied to nonuniform image samples to generate a visuaSy indistir>guishable reconstructed image. 



A warp function relates the {x, y) points in the input 
(original) image to their (V.y) counterparts in the output 
(warped) image. We can define the mapping in terms of 
two funcUons: = X(x,y) and y = Y{x.y). That is, for 
every point {x, in the input image, functions Xand Pre- 
late its correspKjndlng V- and y -coordinates. resp>ectlvely. 
Consider feature p>oints ix^.y^) in the input image and their 
displaced positions (J^.y) in the output image. The MBA 
algorithm can be used to determine a warp function from 
the feature point displacements by constructing two 
smooth surfaces for functions X and Y. One surface for X is 
required to pass through points (x^.y^*^)' while the other 
surface for Y passes through {x^ - • J^) • 

An image warping example is shown Fig. 8. The origi- 
nal image is given in Fig. 8a. Fig. 8b shows the user- 
specified features consisting of points placed on charac- 
teristic parts of the face, such as the profile, eyes, nose, 
and mouth. These feature points are moved to new posi- 
tions to warp the original image. Fig. 8c shows the new 
positions and the deformation due to the warp function 
derived by the MBA algorithm. The warped image pro- 
duced by resampling the original image using the derived 
warp function is given in Fig. 8d. 

The dimensions of the original image in Fig. 8a is 360 x 
243. The number of feature points in Fig. 8b is 242. Two 
functions for the y- and y -components of the warp 
function are derived simultaneously by using the 2D val- 
ues of control points. In the MBA algorithm, we used 
itiq - 2, /\) = 1 and = 360 , = 243 for the size of the 
coarsest and finest control lattices, respectively. Fig. 8 



demonstrates that the MBA algorithm can adequately 
generate smooth warp functions from a few selected fea- 
ture points. 

6.3 Object Reconstruction 

Object reconstruction from a set of scattered 3D points is 
Important in geometric modeling. The MBA algorithm can 
be used to solve this problem when a mapping between the 
scattered fxiints and a 2D parametric space Is defined. In 
that case, each control point in the algorithm is associated 
with a 3D value corresponding to the x-, y-, and z- 
coordinates of a point on the object. The reconstructed ob- 
ject is a 3D uniform cubic B-spllne surface that passes 
through the scattered points. 

Fig. 9 illustrates an example. Fig. 9a shows a panoramic 
image of cylindrical range data of a head, as measured by a 
range finder (Cyberware 4020/PS 3-D Digitizer). The pixel 
intensities depict depth values. We convert the image to a 
3D object by mapping each pixel to a 3D point, as shown in 
Fig. 9b. To derive a set of scattered points on the object, we 
first apply the Sobel operator [24] to the Image in Fig. 9a 
and threshold the result. The 3D points corresp>onding to 
the resulting pixels consist of characteristic parts of the 
object that have high curvature. We also sample a coarse 
regular grid (e.g., one out of every ten points along each 
direction) on the object to handle low curvature regions. 
Fig. 9c shows the set of sampled points. The MBA algo- 
rithm is applied to the samples to reconstruct the object in 
Fig. 9d. 
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TABLE 1 

Performance Data (Time Measurements in Seconds) 





Basic 
MBA 


MBA 


Basic adaptive 
MBA 


Adaptive 
MBA 


Multignd 
relaxation 


Surface construction 


4.18 


0.92 






26.62 


Image reconstruction 


3.86 


2.23 






4.15 


Image warping 


2.14 


1.65 


1.76 


0.98 


7.90 


Object reconstruction 


6.54 


3.77 


6.24 


4.38 


12.23 



The 512 X 170 image in Fig. 9a was used as the 2D 
parametric space to apply the MBA algorithm to the sample 
points. The number of samples In Fig. 9c is 4,112. In the 
MBA algorithm, we used n\, = 2, tIq = I and = 512, 
iij^ = 170 for the size of the coarsest and finest control lat- 
tices, respectively. Fig. 9 demonstrates the potential of the 
MBA algorithm in reconstructing an object from a set of 
scattered 3D points. 

7 Discussion 

In this sectioa we first present experimental results for the 
performance and reconstruction accuracy of the MBA algo- 
rithm. We provide insight Into the role of B-spUne refine- 
ment and adaptive control lattice hierarchy in the measured 
performance. The results are compared to that of multigrid 
relaxation for thin plate splines, an efficient method for 
scattered data interpolation. In all cases, the MBA algo- 
rithm is faster than multigrid relaxation. We then consider 
the use of an affine approximation as a preprocessing of 
multilevel B-spline approximation. Finally, we compare the 
MBA algorithm to hierarchical B-spline refinement, a related 
technique presented for object modeling. 

7.1 Performance 

The performance of the MBA algorithm for the examples in 
this f>a|>er are given in Table 1. Four versions of the algo- 
rithm were tested: 

• basic MBA algorithm (vyithout B-spline refinement) 

• MBA algorithm (with B-spline refinement) 

• basic adapdve MBA edgorithm (without B-splirie re- 
finement) 

• adaptive MBA algorithm (with B-spline refinement) 

The first row was derived from the surface construction 
example shown in Fig. 4b. The other rows were obtained 
from the examples in Section 6. The performance numbers 
in Table 1 were measured In seconds on a Sun SPARC! 0. 

In the basic adaptive MBA algorithm applied to a control 
lattice hierarchy, the adaptive BA algorithm is used to de- 
rive lattices <t>4 for which 16p < |4>4| , where p is the number 
of data points. The final approximation function / is com- 
puted by adding the functions from each control lattice 
<I>4 as In the basic MBA algorithm. When is available as 
a lattice, we use the forward difference method [14], [51] to 
compute . If <I>4 is stored by a set of control points, how- 
ever, ^ must only be computed in the neighborhoods of 
those control points. In that case, we use a lookup table for 
B-spline basis functions to efficiendy compute in those 
neighborhoods. The basic adaptive and adaptive MBA al- 
gorithms differ only in the use of B-spline refinement for 
reducing the coarser control lattices into one equivalent 



lattice. Note that B-spline refinement can be applied only to 
the coarser control lattices that violate the condition 
16p < |<I>4| . This realizes computational savings in evaluat- 
ing the final approximation function. 

The first row in Table 1 shows that the performance of 
the MBA algorithm Is far superior to that of the basic MBA 
algorithm In constructing the surface In Fig. 4b. This dem- 
onstrates the Important role of B-spline refinement In real- 
izing performance gains. In this case, the adaptive MBA 
algorithms were not applied because the data points are 
sparse and can be interpolated with a rather coarse 64 x 64 
control lattice. The second row in Table 1 shows that similar 
performance gains due to B-spline refinement were realized 
in the image reconstruction example. The adaptive MBA 
algorithms were also not applied in this case because the 
large number of data points violated the condition 
16p < 1<I>4| for every control lattice 0^. 

In the image warping example, the adaptive MBA algo- 
rithm was the fastest among all tested. In that algorithm, 
the approximation functions are computed efficiendy on 
the coarser control lattices by using B-spllne refinement. On 
the finer control lattices, the adaptive BA algorithm is used 
and the approxinnation functions are efficiendy computed 
by applying a lookup table for B-spline basis functions to 
the set of necessary control points. In this case, the MBA 
algorithm Is slower because it must apply B-spline refine- 
ment to a large number of control points on the finer lat- 
tices. In general, if the number of data points is small and 
the control lattice size Is large, B-spline function evaluation 
using a lookup table can be more efficient than B-spline 
refinement Hence, the adaptive MBA algorithm outper- 
forms the MBA algorithm on the finer control lattices In 
spite of the overhead of sorting the arrays and to 
compute the necessary control point values. 

The basic adaptive MBA algorithm is slightiy slower 
than the MBA algorithm in the image warping example. 
Although the basic adaptive MBA algorithm efficiendy 
computes the approximation functions at the finer levels, 
the performance loss is due to the computation required to 
densely evaluate the functions at the coarser levels rather 
than applying B-spline refinement The difference between 
the numbers for the two adaptive MBA algorithms shows 
that B-spllne refinement is very helpful for the coarser levels. 

The last row in Table 1 shows that the adaptive MBA al- 
gorithm Is slower than the MBA algorithm In the object 
reconstruction example. In this case, the arrays and o)^ 
are very large because there are a large number of data 
points. Sorting these large arrays requires a computation 
overhead that undermines the efficiency of the adaptive 
MBA algorithm at the finer levels. Actually, the evaluation 
of 2m approximation function using a lookup table at a finer 
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level is not very efficient in this case due to the large num- 
ber of data points. 

In the last column of Table 1. the performance data is de- 
rived by using multlgrld relaxation to obtain thin plate 
splines for the examples in this paper. The muldgrld relaxa- 
tion method has been considered an efficient numerical tech- 
nique to solve a linear system [5]. The method was used in 
early vision problems -to accelerate the^ numerlcai soliition of 
thin plate surfaces that interpolate scattered data points [48], 
[49]. An alternate api>roach using hierarchical basis functions 
or wavelets has been proposed to provide an efficient solu- 
tion to scattered data interpolation with thin plate splines 
[461. [25]. In the approach, a set of hierarchical basis functions 
is used to accelerate the convergence of an iterative numeri- 
cal technique such as corijugate gradient descent 

Let an interpolation function be evaluated on a grid of 
size Mx N. The time complexity of the multlgrld relaxation 
method is 0(M/V). while the approach using hierEirchlcal 
basis functions runs in 0(MN log M/V). When the grid is 
large, we can expect that multlgrld relaxation derives a thin 
plate surface faster than hierarchical basis functions. Hence, 
multlgrld relaxation for thin plate splines was chosen to be 
compared to the MBA algorithm in Table 1. To implement 
the multlgrld relaxation method, we followed the pseudo- 
code in [5]. where the parameters Vq, Vj. and V2 determine 
the number of relaxations. A summary of this approach to 
derive thin plate splines can be fotmd in [321- We used 
V(, = Vj = V2 = 1 for the surface construction example and 
Vq = 1 and V| = V2 = 3 for the other examples. 

The data in Table 1 shows that multigrid relaxation for 
thin plate splines is slower than any version of the MBA 
algorithm. The performarKe of multigrid relaxation 
strongly depends on the size of the underlying grid and 
remains nearly constant regardless of the number of data 
pKtlnts. In the surface construction example, the approxi- 
mation function is evaluated on a 1,024 x 1.024 grid and 
multigrid relaxation is several times slower than the MBA 
algorithm. In the image reconstruction example, the per- 
formance gap is smaller than In the other examples due to 
the large number of data points. 



7.2 Reconstruction Accuracy 

To demonstrate the accuracy of reconstruction by the MBA 
algorithm, we performed experiments with several test 
functions. Given a test fuiKtion /(x. >). we first sampled 
data points from it and applied the MBA algorithm to ob- 
tain an approximation function g. The difference between g 
and f is then measured by computing the root mean square 
error between the function values on a dense grid. 

We selected five out of the six test functions used in 
Nielsons experiment [41) for trlvariate scattered data inter- 
polation. We simpllRed the selected functions to make them 
bivarlate functions. The resulting test functions are 

f,(x, y) = a75ex J --^ 4 ' + 0-75 expi 55 ^^--^ 

4(x, y) = (tanh(9 - 9x - 9y) + 1)/9 . 

f^{x, y) = (L25 + cos(5.4y))/(6 + 6(3jr - l)^) . 

y) = exp[-^((if- 0.5)* + {y-O.sf)^^ , 
f^{x, y) = ^64-8 l((x - 0.5)' + (y - ttS)' )/9 - 0.5 

In addition, we used a bicubic B-spline surface 4 as a test 
function. Fig. 10 shows the test functions. The domain of 
the functions is {(x.y) | 0 < x S 1. 0 < y ^ l}. The approxi- 
mate ranges of the functions are /;: [0.0.1.3], fj: [0.0.0.22], 
^:[0.0.0.38], i^: [0.0. 0.34], l^: [0.0,0.39]. 4:[-0.1. 1.O]. 

Fa- each test function, we used four data sets. MIOO. M500. 
L160. and C160. MIOO and M500 are small and large data sets, 
which consist of 100 and 500 points, respectively. We uni- 
formly sampled 7x7 and 15 x 15 data points for MIOO and 
M500, respectively, while the other points were randomly 
sampled. LIGO consists of data points sampled from several 
lines, which is similar to sampling from the edges of an ob- 
ject. We used eight lines, and 20 points were sampled and 
perturt)ed on each of the lines. In CI 60, the sampled points 
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were clustered, which is analogous to real-world sampling 
[41]. C160 consists of eight dusters, each of which has 20 
points. Fig. 1 1 shows the sampling positions for the data sets. 



(a) M100 



(b) M500 



(c) L160 (d) C160 

Fig. 1 1 . Sampling positions for test data sets; 

For each test function f, we derived an approximation 
function g by applying the MBA algorithm to each of the 
data sets, MIOO. M500. L160, and C160. The root mean 
square (RMS) error between f and g was measured on a 
dense grid over the domain. That is. 



RMS = 



(M+l)(iV+l) 



where = i/M.Yj = J/N, and M = N = 50. To normalize 
the error values, we divided the RMS error by the differ- 
ence of the maximum and minimum values of / over the 
domain. Table 2 shows the results. 

TABLE 2 

Normalized RMS Errors Between Test FuNCTioh4S and 
Their Approximations 





f^ 


^2 


'3 








M100 


0.016 


0.025 


0.013 


0.006 


0.027 


0,016 


M500 


0.001 


0.005 


0.003 


0.0008 


0.007 


0.002 


L160 


0.031 


0.032 


0.042 


0.008 


0.049 


0.022 


C160 


0.082 


0.097 


0.130 


0.086 


0.080 


0.081 



Table 2 demonstrates that the MBA algorithm recon- 
structs a test function quite accurately regardless of the type 
of the function. The normalized RMS error is less than 5 
percent of the range of the function values, except for the 
case of cluster sampling. It is clear from the examples of 
MIOO and M500 that the accuracy Increases when the num- 
ber of data points gets larger. 



7.3 Initial Linear Approximation 

Only at the coarsest level k ■■ 
original data points P in the MBA algorithm. At each suc- 
cessively finer level it > 0, data contains the residual er- 
ror between the intermediate approximation function g^.j 
and the original data P. The BA algorithm applied at level k 
attempts to remove the residual error by deriving a B-spllne 
approximation of that may be added to g^., to yield a 
more accurate approximation . In that case, if the prox- 
imity data set of a control point Is empty, has no ef- 
fect In reducing the residual errors. Then, the most natural 
choice for the value of is zero so that the approximation 
from the coarser levels remains unaffected. When ^„ has a 
nonempty proximity data set Pij is determined by a 
weighted average of the values 0^ computed by (4). Among 
the solutions of (2). we chose (4) to minimize the deviation 
of a B-spllne approximation from zero. This serves to retain 
the approximation from the coarser levels as much as p>os- 
sible while properly reducing the residuals. 

Regardless of the proximity data set of a control point, 
the BA algorithm determines its value so as to minimize the 
changes to the approximation function obtained from the 
coarser levels. This Is Intuitive for all but the coarsest level, 
which must provide a good initial approximation of the data. 
In general, the BA algorithm applied at the coarsest level 
generates a reasonable global approximation of data points. 
In the event that the coarsest control lattice is very fine, many 
of Its control points may have empty proximity data sets and 
be assigned zero. The final app>roximation will thereby tend 
to contain local peaks near the data pcdnts, which may be 
unsatisfactory for many applications (see Fig. 4c). In this 
case, we may improve matters by using a least squares lin- 
ear fit to obtain a better initial approximation. 

We now review how to compute a linear approximation 
function f| from given data points, P = {(Aj.,^^, zj 
|c=l p}.Let 

f_^{x.y) = a^x-^ 32^+ 33. (6) 

where a, . 33 , and % are parameters. We may fit plane f , to 
data P by solving for x in the set of linear equations 
Ax«b. 

where 



A = 


^^1 Yx 
X2 Yz 


\^ 
1 


X = 




andb = 






y? 


i 

> 











We derive the least-squared solution for f^, that minimizes 
the error (z^ - f , {x^^Yc)^^ • solution can be computed 
by simple linear algebra [42] using the pseudoinverse of A: 



X =r (A^A) ' A^b, 



(7) 



where A^ is the transpose of matrix A. The product A^A 
yields a 3 x 3 matrix whose inverse can be computed effi- 
clentiy. The resulting function fj Is the best linear ap- 
proximation to data P in terms of the squared approxima- 
tion errors. 
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(a) given data (b) initial linear approximatibn 



(c) - - 16.m^ = = 64 (d) m„ = = 1.m^ = ^4 

Fig. 12. Reconstructed surfaces using the MBA algorithm with an initial lirtear approximatkm. 



When we use an initial linear approximation, we first 
derive fuhcaon from the original data Fby (7). Then, the 
MBA algorithm is applied to the deviations Pq of f | from P. 
where = {(x,,y,,^°z^) j A°z, - z,- f_i{x.yi}. The final 
approximation function is the sum of ft and the function 
obtained by the MBA algorithm from Pq. With this proce- 
dure» every level in the MBA algorithm, including the 
coarsest level, approximates the deviations between the 
current intermediate function and the data points. This jus- 
tifies the use of a minimal norm solution (4) and assigning 
zero to control points having empty proximity data sets in 
the BA algorithm. 

Fig. 12 demonstrates the role of a least squares linear fit 
as an initial approximation in the MBA algorithm. Fig. 1 2a 
shows the same data given in Fig. 4a. Fig. 12b depicts the 
linear approximation function f j fitted to the data. In Fig. 
12c. the Initialization is applied to the case in Fig. 4c. The 
approximation function In Fig. 12c continues to have peaks 
and valleys near the data points, as in Fig. 4c. Note, how- 
ever, that the function values are no longer zero away from 
the data. They are now derived from the underlying plane 
f 1 . Fig. 12d shows the result of applying the linear initiali- 
zation to the case in Fig. 4b. The resulting approximation 
function remains virtually the same as that in Fig. 4b. In 
this case, the coarsest level in the MBA algorithm generates 
a reasonable global approximation of data p>olnts. and the 
linear fit is not warranted. 

7.4 Comparison to Hierarchical B-Spllne Refinement 

Forsey and Bartels proposed hierarchical B-spllne refinement 
in [161 to gain finer control In editing a surface. The pur- 
pose of the control lattice hierarchy Is to restrict refinement 
to the locality at which ian editing effect is desired. That pa- 
per, however, did not describe how to manipulate the hier- 
archy to interpolate a set of data points. To address this 
problem, they presented a technique in [17], [18] that augments 



B-spllne approximation with hla-archical B-spUne refinement to 
Interpolate a grid of data using a control lattice hioarchy. Their 
method, however, cannot handle scattered data. 

In this paper, we presented multilevel B-spllnes to han- 
dle scattered data. Instead of locally refining control lat- 
tices, we maintain all control points In the hierarchy to ap- 
ply B-spline refinement. This allows us to reduce a se- 
quence of B-spllne functions, defined by a control lattice ~ 
hierarchy, as one equivalent B-spllne function, defined by 
one control lattice. In addition to the performance gains in 
evaluating the app-oximation function, this reduction 
makes the final approximation function suitable for any 
conventional modeling and rendering system. In contrast, 
the hierarchical B-spline refinement technique needs a 
complicated data structure to maintain the local refine- 
ments of control lattices. In this paper, even when an adap- 
tive control lattice hierarchy is used, a control lattice can be 
represented by a simple set of necessary control points. 

8 Conclusion and Future Work 

This paper has presented a fast approximation and inter- 
polation algorithm for scattered miiltivariate data. The al- 
gorithm is based on B-spIine approximation. A B-spllne 
control lattice Is efficiently determined by minimizing a 
local approximation error for each control point. The re- 
sulting -continuous function passes near densely distrib- 
uted data points and interpolates Isolated points. However, 
a tradeoff exists between the shape smoothness and ap- 
proximation accuracy of the function, depending on the 
control lattice density. 

Multilevel B-spllne approximation was introduced to 
circumvent this tradeoff. The algorithm makes use of a hi- 
erarchy of control lattices to generate a sequence of func- 
tions whose sum approaches the desired approximation 
function. In the sequence, a function from the coarsest con- 
trol lattice provides an Initial estimate, which is further re- 
fined in accuracy by functions derived at finer levels. Large 
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performance gains are realized by applying B-spline re- 
finement to reduce the sum of these functions into one 
equivalent B-spline function. 

We also presented a sufficient condition for multilevel B- 
spline approximation to generate an interpolation function 
through scattered data points. Interpolation Is achieved 
when the control point spacing in the finest lattice becomes 
sufTiciendy small relative to the data distribution. We pro- 
posed an adaptive representation of the control lattice hier- 
archy to minimize the memory overhead that may be In- 
troduced at finer levels. 

In his work comparing several scattered data interpola- 
tion methods, Franke used accuracy, visual aspects, sen- 
sitivity to parameters, timing, storage requirements, 
and ease of Implementation as the tested characteris- 
tics of the methods [20]. In Section 7.1 and Section 7.2, we 
showed that the MBA algorithm has good performance 
with respect to the timing and accuracy characteristics, re- 
spectively. In terms of visual aspects, the algorithm gener- 
ates an approximation function with a smooth shape. The 
parameters of the algorithm are the sizes of the coarsest 
and finest control lattices in the control lattice hierarchy. 
The effects of those parameters are quite intuitive as illus- 
trated in Fig. 4. The storage requirement of the MBA algo- 
rithm relates to the size of the finest control lattice, which 
can be controlled by the user. As demonstrated by the 
pseudocodes in this paper, the MBA algorithm is very easy 
to implement. 

The multilevel B-spline approximation can reconstruct a 
planar surface when we apply an initial linear approxima- 
tion to the data, as described in Section 7.3. When data 
points are sampled from a nonplanar surface (e.g., a piece- 
wise bicubic surface), the algorithm does not exacdy recon- 
struct the original surface. However, the algorithm is guar- 
anteed to generate a -continuous piecewise bicubic sur- 
face interpolating the data points, which is expected to be a 
nice approximation of the original surface. On the other 
hand, the algorithm is not afflne-invariant Since data 
points are projected to underlying lattices of varying resolu- 
tion, the relative positioning between the data points and 
the lattices affects the approximation function. In future 
work, we will conduct a complete analysis of the recon- 
struction class of the algorithm and investigate modifica- 
tions to satisfy the afflne-invariance property. 

The scattered data interpolation technique introduced 
here can be applied to other areas of computer graphics. We 
have demonstrated its use for warp generation in image 
warping. This has direct impact on image morphing as well. 
Furthermore, the algorithm may play a significant role In 
data compression. We have demonstrated that high fidelity 
image and object reconstruction is possible from a selected set 
of sparse and irregular samples. Future work will address 
techniques to uniquely determine a minimal set of samples 
necessary to achieve high quality reconstruction within a 
specified error tolerance. 
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